An Improved Forecast of Patchy Reionization Reconstruction with CMB 
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Inhomogeneous reionization gives rise to angular fluctuations in the Cosmic Microwave Back- 
ground (CMB) optical depth r(n) to the last scattering surface, correlating different spherical har- 
monic modes and imprinting characteristic non-Gaussianity on CMB maps. Recently the minimum 
variance quadratic estimator f (n) has been derived using this mode-coupling signal, and found that 
the optical depth fluctuations could be detected with (S/N) 2 ~ fOO in futuristic experiments like 
CMBPol. We first demonstrate that the non-Gaussian signal from gravitational lensing of CMB 
is the dominant source of contamination for reconstructing inhomogeneous reionization signals, 
even with 98% of its contribution removed by delensing. We then construct unbiased estimators 
that simultaneously reconstruct inhomogeneous reionization signals r(n) and gravitational lensing 
potential <f>(n). We apply our new unbiased estimators to future CMB experiment to assess the 
detectability of inhomogeneous reionization signals. With more physically motivated simulations of 
inhomogenuous reionizations that predict an order of magnitude smaller C[ T than previous stud- 
ies, we show that a CMBPol-like experiment could achieve a marginal detection of inhomogeneous 
reionization, (S/N) 2 ~ 0(1) with this quadratic estimator to ~ 0(10) with the analogous maximum 
likelihood estimator. 

PACS numbers: 



I. INTRODUCTION 

Reionization marks the time in which the vast major- 
ity of the hydrogen in the Universe was ionized. When 
and how this process occurred is at present poorly con- 
strained. Current data show that it must have ended 
by z ~ 6 because at lower redshifts there was significant 
transmission in the Lya forest [l| • In addition, the large- 
scale polarization anisotropics in the cosmic microwave 
background (CMB) constrain the mean redshift of reion- 
ization to have been z — 10.6 ± 1.2 @. 

It is believed that the first galaxies in the Universe pro- 
duced the ionizing photons that ultimately ionized the in- 
tcrgalactic gas (e.g. [1]). The morphology of reionization 
and its duration depended on the nature, abundance, and 
clustering of the ionizing sources [1, Q . There are several 
established ideas for how to better constrain the mor- 
phology and the redshifts over which it occurred. These 
include detecting the reionization-induced suppression 
and spatial modulation in the statistics of high-redshift 
Lyman-a emitting galaxies 043 > studying H i Lyman- 
a damping wing absorption from the neutral gas dur- 
ing reionization in the afterglow spectra of high-redshift 
gamma ray bursts [l,[!|[l(|, and directly observing 21 cm 
emission from z > 6 neutral hyd rog en in the intergalactic 
medium (e.g., Furlanetto et al. [1JJ). This study concen- 
trates on using a new technique, first proposed in Dvorkin 
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and Smith [12j, that exploits the non-Gaussianities in the 
CMB sourced by reionization to study this process. 

Inhomogeneous reionization produces several sec- 
ondary anisotropics in the CMB. First, extra tempera- 
ture (and, to a lesser extent, polarization) anisotropics 
are generated from peculiar motion of ionized regions 
during the entire reionization process (l3l - fl9l |. These 
anisotropics are termed the kinetic Sunyaev-Zeldovich 
effect. Second, ionized bubbles scatter the local CMB 
temperature quadrupole, generating fluctuations in the 
polarization at large scales [l5l [20T ] . Finally, the patchy 
nature of reionization would have resulted in the Thom- 
son scattering optical depth to recombination, r(n), de- 
pending on direction [l5[, [IH, [2lT[23| . Such optical depth 
fluctuations act as a modulation effect on CMB fields by 
suppressing the primordial anisotropics with a factor of 
e -T(n)^ correlating different spherical harmonics. Infor- 
mation contained in r(n) fluctuations could potentially 
probe the duration of hydrogen reionization and the size 
of the ionized regions. 

It is well known that gravitational lensing also imprints 
a non-Gaussian signature on the CMB. Minimum vari- 
ance quadratic estimator has been introduced by using 
the coupled modes to reconstruct the projected lensing 
potential [2i - [27j . Recently, Dvorkin and Smith [2 fol- 
lowed similar technique and derived the minimum vari- 
ance quadratic estimator for r(n). Utilizing a toy model 
for reionization, they estimated that the patchy reion- 
ization signal could be detected with (S/N) 2 ~ 100 for 
a CMBPol-like experiment, with beam full-width half- 
maximum (FWHM) of Ofwhm = 4', and noise sensitiv- 
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TABLE I: Description of various reionization models, 
tances are in comoving units. 



Dis- 





Cr 


d(/dz 


Zmfp (M P C) 


T 


Lbox (Mpc) 


N g 


A 


10 





OO 


0.063 


200 


256 


B 


10 


18 


00 


0.112 


200 


256 


C 


30 





10 


0.090 


200 


256 


D 


20 


16 


10 


0.115 


200 


256 



ity At = 1/^k-arcmin. In this paper, we quantify the 
impact of lcnsing induced non-Gaussianities on the re- 
construction of r(n). We show that lensing biases the 
reconstruction of r(n), and as a solution we construct an 
unbiased estimator for r(n) in the presence of lensing. 

The structure of this paper is as follows. Section |TT] 
provides simple estimates for the size of r(n) fluctua- 
tions, and it describes the cosmological reionization cal- 
culations used here to produce r(n) maps. Scction lLUl dc- 
rives the minimum variance quadratic estimator for f (n) 
in the flat sky approximation, and it quantifies the effect 
on lensing on the estimator. In section fVTl we summarize 
our results and discuss the implications. In Appendix |B1 
we describe in more detail our simulations to reconstruct 
r(n) in the presence of lensing. 



II. INHOMOGENEOUS REIONIZATION 

Patchy reionization produced a line-of-sight dependent 
optical depth that can be written as 



r(n) 



a dz 



ot n e (z) [1 + S b (h, z) + 5 x {n, z)] , (1) 



where ot is the Thompson scattering optical depth, n e (z) 
is the average free-electron proper number density, and 
Si, and S x are the over-densities in baryons and in the 
ionized fraction, xi. 

The angular power spectrum CJ T in the flat sky ap- 
proximation can be related to the 3D ionization and den- 
sity field using the Limber approximation: 



dV 2 2 - I \2 

— a <T T n e {z) 



Pxx{^: ) 



+ 2P xS (z,-)+P S s{z,-) 
V V . 



(2) 



where rj is the conformal distance from the observer, and 
Pxy is the cross power spectrum of the over-density in 
X with the overdensity in Y. This power spectrum is 
weighted heavily to the highest redshifts where there was 
reionized gas. The kinetic Sunyaev-Zeldovich effect (kSZ) 
signal from reionization is predicted to be comparable 
to the kSZ signal after reionization. However, the kSZ 
weights by an additional v 2 factor which goes roughly as 
the scale factor, result in its kernel peaking at lower red- 
shifts [lH . This results in a large fraction of the kSZ com- 
ing from after reionization, we can safely neglect the low 



rcdshift part, whereas we expect most of C[ T to originate 
from during reionization. It is also clear from equation 
(|2J) that C[ T from reionization increases approximately 
linearly with the duration of reionization for fixed mean 
redshift of reionization. 

It is likely that reionization occurred in a patchy man- 
ner, with some regions being ionized early on in this pro- 
cess and others remaining neutral until the end, and with 
little gas at intermediate ionization states. This patch- 
iness likely resulted in the ionization fluctuations domi- 
nating over other sources of fluctuation (i.e., P xx 3> P$$ 
on arcmin and larger scales [|[ ) . Even without any knowl- 
edge of P xx other than that reionization was patchy, there 
is an integral constraint on C[ T because if the ionization 
field is zeros and ones / k 2 dk / (2-k 2 ) P xx = a;" 1 — 1, where 
Xi is the ionized fraction. Thus, fixing the reionization 
history and in the Limber approximation, j l 2 dl C'[ T is 
just a single number independent of morphology. This 
constraint shows that the larger the H n regions during 
reionization, the larger the fluctuations in r(n) 1 . 

To estimate CT T , we compute Monte-Carlo realizations 
of reionization in two hundred comoving Mpc data cubes 
using the method developed in [28| for assigning the 
ionization state to boxes with realizations of the linear- 
theory cosmological density field. This method is based 
on the semi-analytic model for reionization in The 
distribution of ionized gas found in the (2&| method is in 
excellent agreement with the results of detailed numeri- 
cal simulations of reionization [28|, |2jj. Thus, we expect 
that the t field from this simulation will be more realistic 
than the analytic model used in the original study of fl2l ] . 
Their model assumed a lognormal distribution of bubbles 
with a distribution that was independent of ionized frac- 
tion. In our calculations, the morphology of the bubbles 
is complicated and their sizes increase dramatically as Xi 
increases. 

The method in [28| that we employ posits that the 
number of galaxies within a region sets its ionization 
state. Namely, a region is ionized if 1 > ( /, where £ 
is a factor that encodes the efficiency that galaxies can 
ionize their surroundings, and / is the total fraction that 
has collapsed into halos with mass > m m j n , where m m j n 
is the minimum halo mass of the sources during reion- 
ization. A point in space is marked as ionized if this 
criterion is met for any smoothing scale centered around 
it (where the smoothing is done with a tophat in Fourier 
space filter). 

Calculating / in detail requires high-resolution iV-body 
simulations to resolve ~ 10 s M Q halos - the smallest ha- 
los that were expected to form multiple generations of 
stars -, while still capturing scales much larger than the 
10 comoving Mpc bubbles. Fortunately, extended Press- 



Since I Ci ~ const., if the peak in the bubble scale is at smaller 
I (i.e. larger bubbles), as l' 2 Ci ~ I , larger bubbles result in a 
higher peak (but at lower / i.e. larger fluctuations). 
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the duration of reionization while the later parameter 
primarily affects its morphology [f| |32[ . 

In particular, bubbles that are larger than mean free 
path of ionizing photons have most of the photons pro- 
duced within them absorbed by dense systems within the 
bubble rather than by diffuse gas at the bubble edge, pre- 
venting further growth [H, |32| . Thus, the parameter Z m f p 
is implemented by setting the maximum smoothing scale 
used to be Z m f p . We generate Monte- Carlo maps for 4 
different reionization models, which are described in Ta- 
ble |TJ All of our models fall within 2a of the best fit 
WMAP r measurement of r = 0.088 ± 0.015 0. 

In Fig.[TJ the top panel shows the optical depth fluctu- 
ation power spectrum of the different reionization models 
described in Table [TJ The corresponding reionization his- 
tory of the four models are also shown in the bottom 
panel of Fig. [JJ Surprisingly, the spectrum of all these 
models is not significantly different: All the models scales 
as 1 2 C[ t /2tt sa constant for 200 < Z < 10000. However, 
the amplitude varies between ~ 10~ 6 - 10~ 7 , owing to 
the different reionization histories. An amplitude of 10 -6 
is still an order of magnitude smaller that the signal con- 
sidered in the previous work of [llj]. It is possible that 
reionization is more extended than in our models. We 
note that the amplitude of CJ T is proportional to the 
duration of reionization. 2 



III. STANDARD QUADRATIC ESTIMATOR OF 
PATCHY REIONIZATION FROM THE CMB 



FIG. 1: Upper panel: Optical depth power spectrum from 
the different reionization models described in Table [U Lower 
panel: Corresponding average reionization history of these 
models. The mean optical depth of each model is consistent 
with the WMAP measurement (r = 0.088 ± 0.015). 



Schechter theory provides a method to calculate / in a 
macroscopic region of size R and overdensity Sr in the 
simulations from just the linear density field [30l l3lj ]. 
Therefore, we can quickly compute the ionization field 
from the linear density field and rather course resolution 
using just Fast Fourier Transforms. Our calculations take 
just minutes on a single CPU for the 256 3 grids used here. 

There is significant uncertainty in the properties of the 
first sources and sinks that were responsible for reion- 
ization. All of our models assume that the ionizing lu- 
minosity is proportional to the collapsed mass in halos 
above 10 8 M© (approximately the minimum mass thresh- 
old where the gas can cool by atomic transitions and 
form dense structures). To explore the allowed param- 
eter space of this process, we model reionization with 3 
parameters: the ionizing efficiency of a halo (7 at 2 = 7, 
its derivative with redshift dC^jdz (assumed to be inde- 
pendent of z), and the mean free path of ionizing photons 
to be absorbed by an over-dense sink Z m f p within an ion- 
ized region. The first two parameters primarily affect 



The observed CMB temperature and polarization 
Stokes parameters in the presence of inhomogeneous 
screening caused by patchy ionizated regions are 

T(n) = e- ST ^f(h), 
(Q±iU)(n) = e-* r W(Q±iU)(n), (3) 

where tildes signify the CMB Stokes parameters for a uni- 
form reionization history with constant factor e~ T spa- 
tially modulating the observed CMB fields. We take r 
as the mean of optical depth and 5r(n) as the line of 
sight dependent optical depth fluctuation field. We work 



2 Recently it was shown that the velocity difference between the 
baryons and dark matter that is imparted up until recombina- 
tion and decays away thereafter, can suppress the formation and 
baryonic accretion of the < 10 6 Mq halos that harbor the first 
stars [33I , l34ll . The standard paradigm is that these halos did not 
contribute significantly to reionization [Till , but they may have 
ionized the intergalactic medium fractionally. Different regions 
in the Universe have different velocity offsets, with the coherence 
length of this difference being hundreds of Mpc. Even if these 
first stars just fractionally ionized the Universe, this large-scale 
modulation of the velocity difference could lead to larger fluctua- 
tions in r(n) (and peaking at I ~ 100) than in the models we have 
considered. Thus, we point out that there remains the possibility 
of generating a larger signal than in the models considered here. 
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TABLE II: Minimum variance filters for optical depth estimator f (n) and lensing potential estimator 4>(n) 
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TB 
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EE 
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EB 
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Via) 


BB 
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in the flat-sky limit where scalar fields such as the CMB 
temperature T and a complex field (<Si ± 1S2 )(n) of spin 
±s can be expanded in the Fourier basis as 



T(l) = / dnT(n)e- 4l - A , (4) 
[«Si±i5a](l) = {±l) s J dh[S 1 (h)±iS 2 {h)]e Ts ^e- ll - & , 

where tpi — cos _1 (n • 1). The complex field (Q ± if/)(n) 
is a spin ±2 field, whose Fourier harmonics are referred 
as (E±iB)(l). 

Since the differential optical depth fluctuation is al- 
ready constrained to be small, we work out the effects 
to first order in <5r(n). We use r(n) rather than 5r(n) 
to specify the fluctuations of optical depth for short. It 
is simple to show that patchy reionization induces mod- 
ulations in observed CMB fields that is proportional to 
r(n) to the first order. The effect of such a modulation is 
to correlate different CMB modes in Fourier space. The 
correlations can be compactly written as 



(X(l 1 )X / (l 2 ))c MB = /ix'(li,l2)r(l) 



(5) 



where X, X' =T,E,B,l = h + 1 2 , f X x>( 1 , l ') is S iven in 
Table HH and (...)cmb signifies an ensemble average over 
CMB realizations with fixed r(n) field. 

The presence of r(n) field breaks the rotational sym- 
metry of the CMB field, correlating different modes 
which are not correlated assuming a Gaussian CMB 
field. Following (24|, we construct a minimum variance 
quadratic estimator fxx'(n) for r(n) field, or fxx'(l) for 
t(1) in Fourier space. 



fxx'O) = N XX '(l) 
where 1 = I2 + li and 
d 2 h 



d 2 h 
(2tt)2 



[X(h)X'(h)]F xx ,(h,h), 



N XX >{1) 



fxx>(h,h)F xx >(h,h) 



(6) 



We derive the optimal Fxx' by minimizing the variance 
of (fxx'(l)fxx' (!'))• For XX' = EE, BB, and TT, 



F X x(h,h) 



For XX' = TB and EB, 



F X x>(h,h) 



fxxQuh) ^ 

e^^~~yX X ,t ^~jX X ,t 
ll h 



fxx>(h,h) 

h h 



where 



qXX ,t ^jX X _j_ ^jXX , 



(7) 



(8) 



(9) 



and Cf X ,n is the noise power spectrum. We assume the 
detector noise is Gaussian and isotropic, to be known a 
priori. Furthermore, we assume a symmetric Gaussian 
instrumental beam so that the noise power spectrum is 



c XX,n = A 2 fe i 2 OL hm /(8 ln2) j 



(10) 



where Ax is the instrument noise for temperature (X = 
T) or polarization {X = E,B), and 0f w hm is the full- 
width half-maximum (FWHM) of the Gaussian beam. 
We assume a fully polarized detector, for which Ae.b = 
V2A T - 

The variance of the minimum variance quadratic esti- 
mator is 

(fxx'(h)Txx'(h)) = (2tt) 2 < 5(1 1 +h){Cr + N XX '(l)}, 

(11) 

where N XX '(J) gives the dominant contribution to the 
variance for the EB and TB estimators. 



IV. LENSING CONTAMINATION IN r(ii) 
RECONSTRUCTION 



The optical depth estimators described in the previous 
section neglect the effect of CMB lensing. In reality, both 
the CMB temperature and polarization fields are gravi- 
tationally lensed by inhomogeneities in the matter distri- 
bution between the last scattering surface and z = 0. In 
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Input Reconstructed Difference 




FIG. 2: Impact of lensing on the reconstruction of the optical depth fluctuation field r(n). The upper panels assume all the 
mode-coupling of CMB maps are generated by modulations of patchy reionization only. The upper left panel shows the input 
r(n) map used to modulate the CMB fields, the upper middle panel shows the reconstructed r(n) map from the CMB fields 
by applying the quadratic EB estimator, and the upper right panel shows the difference between the input /reconstructed r(n) 
maps. The lower panels show the same quantities as the upper panels but with the lensing effect on the CMB maps included. 
As is clear from the lower middle and lower right panels, additional mode-coupling of CMB fields due to lensing contaminates 
the T~(n) reconstruction. In this plot, we consider the lensing signal to be only 2% of the fiducial value (i.e. the deflection 
angle power spectrum used here is Cf d /50) to approximate the residual lensing signal after applying delensing procedure on the 
observed CMB maps. The reconstructed maps were averaged over 1000 CMB realizations for a fixed optical depth fluctuation 
r(n) field. Each map is 6x6 square degrees. 



this section, we show that lensing significantly bias the 
r(n) reconstruction. 

Both the r(n) field and the projected lensing potential 
(j> len (n) can generate non-Gaussianity by mixing modes 
and break the rotational invariance. This effect can be 
detected statistically by searching for the characteristic 
four point correlations. If the f(n) estimator derived 
in the previous section were applied to the lensed CMB 
maps, it would also pick up significant spurious signal 
produced by lensing. 

We now quantitatively calculate the lensing bias to 
the f(n) estimator. Lensing simply deflects the path of 
CMB photons from the last scattering surface resulting in 
a remapping of the CMB temperature/polarization pat- 



tern on the sky. The deflection angle d(n) is related to 
cj) len (h), the lensing gravitational potential as 

d(n) = V0 /e "(n). (12) 

The lensing potential (j> len (h) is given by 

ie "(A) = -2 f°dr ^\°~ r) r $(r,n), (13) 
J Q dA{r)dA{r ) 

where cLa is the co-moving distance along the line of sight; 
ro is the comoving distance to the surface of last scatter- 
ing, and $ is gravitational potential [27j ■ 

Similar to the effect of screening from patchy reioniza- 
tion, a lensing potential mode with wavevector 1 mixes 
the two polarization modes of wavevectors li and I2 = 
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FIG. 3: Comparison of the biased reconstruction of patchy 
reionization power spectrum C[ T with the input fiducial 
patchy reionization induced power spectrum CJ T . The solid 
red line shows the CJ T of the reionization scenario we chose. 
The dot magenta curve shows the biased estimate for CJ T , 
which was calculated using Eq. Q5) . The lensing bias is neg- 
ative at low I but positive at high I. The minimum of the bias 
at / ~ 200 is where the sign changes. The green dashed curve 
is the same quantity extracted from simulations after averag- 
ing over 1000 realizations. Our analytic expression matches 
well with the simulation. The blue short-dashed curve shows 
the reconstructed CJ T after delensing is applied. The de- 
lensing was performed using the quadratic minimum variance 
estimator for projected lensing potential. 



1 — li . Taking the ensemble average of the CMB fields for 
the fixed </>(n) field, similar to Eq. ((5J), one gets 



(X(1 1 )X'(1 2 )) C MB = W e "(l) 



(14) 

The form of filters fxx'O) f° r different combinations of 
CMB field X and X' are given in TableHH The major dif- 
ference between the filters of lensing potential estimator 
</>(ri) and those for the f(n) estimator is the additional 
factors of ~ I 2 that appear in those for lensing owing 
to the differential in Eq. (|T2|) . This differential nature 
of lensing significantly suppress the level of lensing esti- 
mator noise since I ^> 1 and iVj^ s , (Z) is approximately 
proportional to l~ 2 (see Eq. ©). 

Let us consider a CMB sky that has been modified 
by both inhomogeneous reionization and lensing. Sup- 
pose we want to reconstruct r(n) assuming that all of 
the non-Gaussianity is from patchy reionization, which 
is equivalent to applying the f(n) estimator with filters 
designed to optimally reconstruct r(n). In this case, the 
estimator measures 

(fxx'( 1 ))cMB = r(l) 
d 2 h 



+ N xx ,(l) 



(27T) 



zfxx'Fxx<{h,hW en {\), (15) 



bias 

where Nxx> is given by Eq. ©. The first term on the 
right hand side is the desired signal, and the second term 



is a bias that owes to lensing. Note that f xx , is the 
lensing filter (see Table HI)) and F xx ,(li,l 2 ) is given by 
Eq. and ©. 

We simulate a patchy reionization induced r(n) field 
(model B in Table H and modulate the CMB fields by 
the r(n) field accordingly. We compare the reconstructed 
r(n) with the input field in Fig. [2] (sec Appcndix[B]for de- 
tails of the simulations). The r(n) estimator is unbiased 
if primordial CMB fields were unlensed and only affected 
by patchy reionization. However, in the presence of lens- 
ing, the reconstructed r(n) deviates significantly from 
the fiducial signal. The lensing induced non-Gaussianity 
is rougly an order of magnitude larger than the patchy 
reionization induced non-Gaussianity. 

The lensing induced non-Gaussianity could be reduced 
by applying lensing estimator 0(n) to reconstruct the 
lensing potential, and then "remap" the observed CMB 
photons given the reconstructed 4>(h) and Eq. (jT2j) . This 
process of subtracting the lensing effect from CMB is re- 
ferred to as "delensing" (see Ref [3f| for a review). To 
investigate the lensing bias after applying this delens- 
ing procedure, we assume the residual lensing potential 
power spectrum is only 2% of the fiducial value. The de- 
lensing fraction taken here is smaller than the predicted 
delensing fraction for future CMB experiment, using lens- 
ing maps either externally reconstructed from large scale 
structure/CMB temperature or from CMB polarization 
itself [36[ . We find that even after delensing on the CMB 
map, the reconstructed r(n) field is still significantly con- 
taminated by the residual lensing signal, as shown in 
Fig. El 

In Fig. [3j we show the reconstruction of optical depth 
fluctuation power spectrum CJ T , and compare with the 
input power spectrum CT T ■ Again we choose Model B 
for the reionization simulations, which has the highest 
level of r(n) fluctuations. We find that the lensing in- 
duced spurious signal dominates over the fiducial signal 
by ~ O(10 - 100), especially for I > 200. The theoreti- 
cal prediction for the spurious patchy reionization signal 
from lensing which is calculated by Eq. (|T51) . matches 
well with C[ T from the simulation. Finally we show that 
even after applying the delensing procedure with lensing 
quadratic estimator (24[, the reconstructed C[ T is still 
biased by a factor of ~ 10. As we show in Fig. [3j the 
lensing induced C[ T has two bumps one peak at large 
scale I ~ 50 and the other peaks at small scale I ~ 1000. 
It is caused by the lensing bias given by Eq. (p~5|) is neg- 
ative at low I and positive at high I with a transition at 
I sa 200. This sign change is because the lensing bias 
involves the product of the lensing and tau filters [see 
Eq. (|15[) ]. The product contains a mode coupling term 
1 • li which is caused by the derivative nature of lensing 
and gives the negative contribution at low I. Physically, 
the lensing of CMB does not generate new power in the 
CMB fluctuations, it only move power from large scale 
to small scales [231 - We note that in principle lensing 
reconstruction is also biased by the patchy reionization 
induced non-Gaussianity, however since lensing signal is 
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much larger than the patchy reionization signal, we don't 
expect a significant comtanimation from patchy reioniza- 
tion to lensing esitmation. 



V. RECONSTRUCTING PATCHY 
REIONIZATION 

A. Unbiased Estimator 

This section constructs an unbiased estimator for r(fi). 
As with the quadratic estimator discussed in the previous 
section, among all the six estimators, the EB estimator 
has the highest S/N ratio, thus we focus on EB estimator 
in this section. For each multipole I we can define a 2- 



by-2 Fisher matrix F, 



ad 



F, 



ad 



(2tt)' 



where a and j3 run over r and </>. The element (i^" 1 ))^ 3 



of the inverse of Fisher matrix gives the variance of C l 
Hence the variance of CT T is: 



a/9 



N(l) = {¥' 



F rrp00 _ (pry 



(16) 



This is the Gaussian noise term of the unbiased f (n) 
estimator, which we will use to calculate (S/N) 2 in Fig. 2] 
Starting from biased estimator f (n) and </>(n) [24|, we 
have 

pj4 



(f(l))cMB - r(l) + -^0(1) 



(0(1)) 



CMB 



F <I>T 

0(1) + ^r(l) 



(17) 



One can then solve above equations for r(l) and </>(!) 



r(l) 
0(1) 



F^FrimhMB - F^Ft* (0(1)) cmb 



(K 4 



F ^ Frm) 



CMB 



F r F i 



(pry 



(18) 



This estimator although unbiased is not a minimum vari- 
ance estimator. In next subsection we compare the vari- 
ance (Gaussian noise) of the minimum variance quadratic 
estimator with the variance of the unbiased estimator and 
show that there is only a marginal increase in the variance 
of the unbiased estimator in comparison to the variance 
of the minimum-variance estimator. 

Maximum Likelihood Estimator: Given that the B- 
mode polarization is well mapped, Hirata and Seljak 
[26j found that for lensing reconstruction the maximum- 
likelihood estimator (which reduces the estimator noise 
from lensing) allows significantly better (S/N) 2 than the 
quadratic estimator. 



Following Rcf. [26j, the lensing maximum- likelihood 
estimator can be generalized to construct a unbiased 
maximum-likelihood estimator for r(n) in the presence of 
lensing. The variance of the maximum-likelihood estima- 
tor for r(n) is the same as that for the quadratic estima- 
tor r(n) with one exception — for maximum-likelihood 
estimator the denominator of Eq. ([7]) and Eq. ([5]) con- 
tains the unlensed CMB power spectrum, whereas the 
quadratic estimator noise contains the lensed CMB power 
spectrum (26|. The estimator noise of r(n) reconstruc- 
tion would no longer be saturated because of the lensed 
CMB power spectrum. Conceptually, the lensing or 
patchy reionization induced B-modes can be iteratively 
cleaned from the map, therefore we are able to reduce the 
post-cleaning B-mode power spectrum and thus reducing 
the noise in the f(n) estimator. Our fundamental ability 
to clean the map is bounded by the sum of the unlensed 
CMB B-modcs and the instrumental noise. 



B. Forecasting the Detectability of Patchy 
Reionization 



The signal-to-noise for the detection of patchy reion- 
ization signal can be written as 



5n* 



%£(2I + 1)| 



N(l) 



(19) 



where f s k y is the sky fraction; C[ T is the fiducial patchy 
reionization power spectrum, and N(l) is the leading or- 
der Gaussian noise of an estimator, given by Eq. (fT6| for 
the unbiased quadratic estimator and given by Eq. © 
for the biased minimum variance quadratic estimator. 

In Fig. |4l dashed-lines (the lower two curves which 
almost overlap) compare the (S/N) 2 of the biased and 
unbiased quadratic estimators. The left panel is for 
the CMBPol like experiment with with noise At = 
-arcmin and beam FWHM Of w hm — 4 arcmin. The 
right panel is for the reference experiment with noise 
At = 0.2/ik-arcmin and beam FWHM 9f w h m = 1 ar- 
cmin. As is clear from figure the (S/N) 2 of unbiased es- 
timator is only slightly lower than the (S/N) 2 of biased 
estimator for both CMBPol like experiment and the refer- 
ence experiment. In another word, the variance of the un- 
biased estimator is only marginally more (percent-level) 
than the variance of the minimum- variance quadratic es- 
timator. The reason for this is easy to understand — the 
contribution to the variance from the spurious r(n) sig- 
nal produced by lensing is much smaller than the intrinsic 
f(n) estimator noise. 

In Fig. [U dotted-lines (the upper two lines which 
almost overlap) compare the (S/N) 2 of the bi- 
ased/unbiased maximum- likelihood f(l) estimators. For 
CMBPol- like experiment, the maximum-likelihood esti- 
mator can get (S/N) 2 about a factor of 10 higher than 
the quadratic estimator. 



8 



in 

CD 
W 

'o 



D 
c_ 
cn 



4 Quadratic Estimator 

U Unbiased Quadratic Estimator 

Maximum Likelihood Estimator 

Unbiased Maximum Likelihood Estimator 



1 o 2 - 



1 o° - 



00 



EB estimator 

FWHM = 4,0 orcrnin 

A T = 1 -O^K-orcmin 



1 000 
Angular multiple 



CO 

CD 

m 

'o 



1 4 - 



1 o 2 - 



D 
cn 



10°- 



EB estimator 
FWHM=1.0 arcmin 
A T =0.2/xK-orcmin 



00 1000 
Angular multiple 



FIG. 4: The cumulative (S/N) 2 for the EB estimator as a function of maximum multipole I. The left panel is for a CMBPol- 
like experiment with a beam of Bfwhm = 4' and the noise sensitivity A p = l^I<^arcmin. The right panel is a more sensitive 
experiment with 8fwhm = 1' and A p = 0.2/iT^-arcmin. In both panels, we use the optical depth power spectrum CJ T of 
model B in Fig.l. The green/red dotted curves (lower two curves) which almost overlap are the (S/N) 2 of the biased/unbiased 
quadratic estimators, whose noise level is given by Eq. © and (JTBJ) respectively. The magenta/blue dashed curves (the top two 
curves that nearly overlap ), are the (S/N) 2 of the biased/unbiased maximum-likelihood estimators. The maximum-likelihood 
estimator is calculated from the minimum variance quadratic estimator except with the lensed CMB power spectrum replaced 
by the primary CMB power-spectrum (without the presence of lensing or patchy reionization effects) as suggested in [2(| . 
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FIG. 5: Dependence of the total (S/N) on instrumental sensitivity, A p , (left panel) and beam size, 0fwhiv 
plot EB unbiased quadratic estimator as an example and calculate cumulative (S/N) 2 up to l ma x = 3000. 
beam size with 1, 2, 4, and 6 arcmin respectively to calculate the (S/N) 2 dependence on instrumental sensitivity. Right panel, 
we fix A p with 0.5, 1.0, 1.5, and 2.0/iK-arcmin to show the beam size dependence of (S/N) 2 . The approximate fitting functions 
are given in Eq. ([20l . 



In Fig. El wc show the total (S/N) 2 from the unbiased 
quadratic estimator as a function of instrumental beam 
size and detector sensitivity respectively. We find that 
the (S/N) for a constant l 2 Cf T / (2ir) and for experiment 
with Ap > 1/xk-arcmin can be approximated as 



N 



;i 2 cr/2n 

^ J sky 



x exp 



5 x 10- 6 

— 0.099fwhm 
V 



cxp 



-0.56A, 



1/j.K — arcmin 



(20) 



The (S/N) 2 is more sensitive to the instrumental sen- 
sitivity rather than the beam size. For a reference 
pathcy reionization signal l 2 CT T /2ir = 5 x 10~ 6 , for a 
CMBPol-like or COrE-like [37| experiment wc expect 
(S/N) 2 ~ 0(1). For the future ground-based experi- 
ments such as the POLAR Array with A p = 1.41/iK- 
arcmin, Ofwhm = l'j and sky coverage 100 deg 2 , we 
expect (S/N) 2 ~ 0(0.01). 
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VI. SUMMARY AND DISCUSSION 

Reionization marks the epoch in which the vast ma- 
jority of the hydrogen in the Universe was ionized since 
cosmological recombination. When and how reionization 
occurred is at present poorly constrained. In addition to 
pinning down the epoch of this cosmic phase transition, 
constraints on the reionization history provides us infor- 
mation about the formation of early galaxies. Inhomoge- 
neous reionization would have generated fluctuations in 
the Thomson scattering optical depth r(n) among differ- 
ent lines of sights at the level ~ 10~ 3 . These modulations 
would modify the primordial CMB temperature and po- 
larization anisotropics by inducing a directionally depen- 
dent screening. Such screening couples different modes of 
CMB, converting i?-modes to -B-modes, and introduces 
non-Gaussian signals. 

In this paper, we used a technique that exploits the 
non-Gaussianities in the CMB sourced by reionization 
to study this process, as first proposed in Dvorkin and 
Smith [121 ]. We have introduced the the minimum vari- 
ance quadratic estimator in an intuitive flat sky limit 
and compared it with the estimator for lensing poten- 
tial reconstruction [24| . Lensing induced non-Gaussian 
features would produce a spurious r(n) signal that is 
at least an order of magnitude higher than our semi- 
analytical models predict from patchy reionization. We 
showed that ignoring the lensing contamination would 
significantly bias the reconstruction of optical depth fluc- 
tuation field r(n). Even after applying a delensing pro- 
cedure that used the minimum variance quadratic esti- 
mator for the lensing potential <j>(n), the residual lensing 
bias on the f (n) estimator was still comparable with the 
fiducial value. As a solution, we constructed an unbi- 
ased estimator to simultaneously reconstruct r(n) and 
the lensing potential <p(n) such that the estimate of r(h) 
is not biased by lensing. We found that the S/N of the 
unbiased estimator is only degraded at the percent level 
compared to the original biased r(n) estimator. 

We studied the detectability of patchy reionization 
by considering more detailed i~(n) fields using semi- 
numerical reionization models, which unfortunately yield 
an order-of-magnitude smaller signal than previously 
considered [l2j|. As a result, we found that with the 
unbiased estimator, a CMBPol-like experiment could 
achieve a marginal detection of patchy reionization with 
(S/N) 2 ~ 1 — 10. We characterized the estimator noise 
for various instrumental properties. We find that the 
S/N is only weakly sensitive to the FWHM of detector 
beam with a factor of ^2 degradation of (S/N) 2 by in- 
crease FWHM from 1' to 6'. While the (S/N) 2 decreases 
by a factor of ^2 by increase instrumental noise from 0.5 
to 2 /ik-arcmin. Similar scaling with instrumental char- 
acteristics have been quantified for lensing reconstruction 

in m. 

Large scale CMB fields are also modulated by smaller 
scale T(n) fluctuations due to patchy reionization. As we 
construct the f(n) estimator in flat sky limit, we ignore 



the patchy reionization signal from large scale E / B-mode 
which is generated via scattering of the local CMB tem- 
perature quadrupole by ionized bubbles. The S/N will 
be increased by a factor of ^2 by considering such signal 
on large scales [L2| . 

Although the predicted S/N for a patchy reionization 
detection is only marginal for a CMBPol-like experiment, 
one can cross-correlate with other cosmological data sets 
that are sensitive to the properties of patchy reionization. 
The same population of ionized bubbles would not only 
induce line of sight dependent optical depth of CMB, 
but also correlate with the distribution of galaxies or the 
redshifted 21cm signal. 

At large scales, it is expected that the distribution 
of galaxies correlates well with the neutral gas distribu- 
tion [28[. One can estimate the (S/N) T for a patchy 
reionization detection as 

(S/N) 2 T = l 2 chm .(C^ aT ) 2 /(Cn 2 (21) 

where Z cna r is the characteristic multipolc that con- 
tributes to the S/N (7 char ~ 10 3 ), and C[ f is the f(n) 
estimator variance. The scaling factor Z^ nar is an esti- 
mate for the number of modes that are contributing to 
the signal (the result really does not rely on the fraction 
of the sky r(n) is estimated). 

An estimate for signal-to-noise that can be obtained in 
cross correlation (S/N) Tg is 

(S/N) 2 

~ ''char /sky /rcion (C c har) / (CJ T C gg ) 

= ^char /sky /rcion T C cnar /C ; 

— ^char /sky /rcion r 2 (S/N) T , (22) 

where / s k y and / re ion are the fraction of the sky and reion- 
ization over which surveys overlap, r is the cross corre- 
lation coefficient of galaxies and the r(h) field over the 
same projected volume as the galaxy survey (r ~ 1 on 
large scale). 

Noting that / s k y < 10 -4 is the current size for z ~ 7 
galaxy surveys, it would take a very ambitious survey 
to enhance the r signal in cross correlation compared to 
in the auto-power. Correlating with the diffuse back- 
ground light from early galaxies - the cosmic infrared 
background - is a related and intriguing possibility since 
then /sky /rcion ~ 1 (although, lower redshift emission 
may be a significant noise source in this case) and may 
deserve further study. 

The final possibility that we discuss is cross correlat- 
ing with a survey of redshifted 21cm emission from in- 
tergalactic neutral hydrogen. Such surveys do span a 
significant fraction of the sky and the first generation 
of such endeavors will be in a noise-dominated reg ime in 
which they could benefit from cross-correlation [38( (Note 
that cross correlating with r would be of little interest if 
there existed high S/N 21cm maps). However, redshifted 
21cm analyses remove the modes with small line-of-sight 
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projected wavevectors in the act of foreground cleaning, 
which unfortunately are the modes that contribute to the 
r signal Thus, there would be little signal in this 

cross correlation. 
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Appendix A: Reconstructing Inhomogeneous 
Reionization r(n): Simulation Pipeline 

Our simulation pipeline of optical depth r(n) field re- 
construction follows the procedure in |40|, and we mod- 
ified the code developed for lensing reconstruction in 
[4ll I42I ]. First, we generate primordial CMB polariza- 
tion <3 p ™(n) and U p "(h) maps as Gaussian realizations 
of CMB power spectrum. We choose a standard fidu- 
cial model with a flat ACDM cosmology, with param- 
eters given by Q b = 0.045, il c = 0.23, H = 70.5, n s = 
0.96, n t = 0.0, and r = 0.08. We calculate the theoretical 
lensed and unlensed CMB power spectrum from publicly 
available code CAMB [Hj]. The primordial CMB polar- 
izations maps are then transformed according to Eq. Q 
to include the effect of patchy reionization. The r(n) field 
was generated from a reionization simulation described in 
Section HU 

To include the effect of lensing, we generate a realiza- 
tion of lensing deflection field d(n) and transform the 
CMB fields Q(n) and U(h) to Q(n) and U(h) according 
to 



(Q±iU){n) = (Q±«E/)(n + d(n)). 



(Al) 



The deflection angle at each point n is calculated by tak- 
ing the gradient of the lensing potential. The lensing po- 
tential power spectrum is generated using CAMB which 
was run with nonlinear corrections using halofit [43| . 

Since we want to quantify the lensing contamination, 
we have several pipelines with different level of lensing 
signal being removed. We define a de-lensing factor a, as 
C deien,<f,<f, =c theory,<f,<f,/ a ^ where a = 1 correspond to no 

de-lensing, a — > 00 corresponds to perfect delensing, we 
use a = 50 for Fig. ©. 

We then Fourier transform CMB polarization maps to 
get E(l) and B(l) maps. Finally we multiply CMB E(l) 
and -B(l) maps by Gaussian beam in Fourier space and 
add instrumental noise. 

More specifically, we closely follow Hu et al [4(| to re- 
write the r(n) estimator which is more efficient to eval- 
uate computationally. We re-write the estimator in real 
space 

ff B = -N EB J d 2 n e --' 1 Rc{[G £;B (n)L B *(n)]} . 

(A2) 

The field G EB is built from the observed £"(1) field (in- 
cluding contributions from lensing and patchy reioniza- 
tion) as 



c EB — 

1 (cr+N EE ) 



E{\)e 



2iip\ 



(A3) 



and L B is given by 



L? = 



B(l) 



{C{ 



BB 



N BB ) 



(A4) 



N EB is a normalization coefficient which ensure the 
unbiansdness of the estimator [4(|. The results of our 
simulations are shown in Fig. (|2|). 



Appendix B: Unbiased Minimum Variance 
Quadratic Estimator for Patchy Reionization 

This section discusses the quadratic estimator for the 
patchy reionization induced optical depth fluctuation 
field r(n) in a more general context, demonstrating that 
the estimator f (n) used in the text is the minimum vari- 
ance estimator in the limit that the signal-to-noise ratio 
in lensing estimator <^(n) is much higher than in f(n). 

In Fourier space, the general unbiased quadratic esti- 
mator for </>i and t\ is 

h = J2 - li)X'(k) - Tr(Qf xx ,) 0, 

U 

= PhX(\ - U)X'(U) - Tr(Pf xx ,) n (Bl) 

u 

where Qii an d Pu are some weighting functions, f xx , 
and f XX i is the same as in Eq. (|5|) and (fT4]) . The 
sum does not include li = and we are using Tr(X) as 
shorthand for J2 U X h- Noting that (X(l - li)X'(k)) = 

fxx' (1 — li 7 h)<fii+ fxx'0-~^i' li) r ij we can write the above 
equation as (if we substitute the unbiased estimator <f>\ 
and f\ on the R.H.S.) 



A x 



n 



WX'(li) 



where 



A=l 1 Tr (Qfxx>) 
MPfxx>) 1 



(B2) 



(B3) 



Thus, the general unbiased quadratic estimator for t\ 
alone is 

n = Y< [[ A_1 ]n Qu + [ A_1 ]i2^i,] r(i - i^Tfo) (B4) 



We want to derive the weighting functions Q\ t and P\ t 
that give us the unbiased minimum variance estimator 
fi. The estimator variance is 

var[f,] = (If.l 2 ) 

= 2j^X h C h X u Oi. h , (B5) 



where 



(B6) 
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To derive the minimum variance estimator, we want 
to minimize var[fi] subject to the conditions that 
Tr[Pf xx ,] = 1 and Ti[Qf xx ,] = 1. Rather than go 
through this exercise, let us note first that at relevant 
multipole f XX i <C fxx' because of the factor of ~ I 2 that 
contributes to f xx , . Let us also note that the weighting 
function Pi i that is optimal for simultaneously estimat- 
ing fa with t\ should be nearly identical to the minimum 
variance quadratic weighting for estimating just fa be- 
cause T\ is a weak contaminant of lensing. Second note 
that KA" 1 ]^ ex Tr[Pf xx ,] < 1 (since Tr[Pf* x ,] = 1), 
and thus [A -1 ] 12 has magnitude that is much less than 
that of [A _1 ] n . Not only [A _1 ]n > [A _1 ] 12 , but note 



the scaling Q\/ P\ ~ I 2 ^> 1, therefore, we are justified in 
ignoring the second term in Eq. (|B6[) and one can show 
that the minimizing Eq. (|B5[) subject to the constraint 

T±[Qfx X >] = 1 y ields 

Q U = nC-'fxx'C^fxx^C^fxx'i^-WCi-X, 

(B7) 

which yields the identical estimator to that used in the 
text. Furthermore, because the variance of fi is domi- 
nated by Qiij this explains why our unbiased estimator 
that accounts for fa yields a result that is not much dif- 
ferent than the biased minimum variance estimator. 



